COVID-19 US

Prepare county level data

Read and format personality data

df_us_pers <- read_csv('timeseries_usa_county_march1_april_09.csv')

df_us_pers <- df_us_pers %>% select(countyfips, open, sci, extra, agree, stabil) %>% 
  mutate(stabil = 6-stabil) %>%
  dplyr::rename(county_fips = countyfips,
         pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = stabil) %>% 
  distinct() %>%
  mutate(county_fips = as.character(county_fips))

df_us_pers

Read and format prevalence data

df_us_prev <- read_csv('USA_timeseries_prep_2005.csv')

df_us_prev <- df_us_prev %>% 
  select(fips, date, rate) %>% 
  mutate(date = as.Date(date, "%d%b%Y")) %>% 
  rename(county_fips = fips, 
         rate_day = rate) %>%
  mutate(county_fips = as.character(county_fips))

df_us_prev

NA

Read and format county level controls

df_us_ctrl <- read.csv('controls_US.csv')

df_us_ctrl <- df_us_ctrl %>% select(-county_name) %>% 
  rename(county_fips = county) %>%
  mutate(county_fips = as.character(county_fips))

df_us_ctrl

NA

Read and format social distancing data FB

fb_files <- list.files('../FB Data/US individual files/Mobility/',
                       '*.csv', full.names = T)

df_us_socdist <- fb_files %>% 
  map(read_csv) %>% bind_rows()

df_us_socdist <- df_us_socdist %>%
  select(-age_bracket, -gender, -baseline_name, -baseline_type, -polygon_name) %>%
  rename(date = ds,
         county_fips = polygon_id,
         socdist_tiles = all_day_bing_tiles_visited_relative_change,
         socdist_single_tile = all_day_ratio_single_tile_users) %>%
  mutate(county_fips = as.character(county_fips))

df_us_socdist$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-03-01" "2020-03-16" "2020-04-01" "2020-04-01" "2020-04-17" "2020-05-03" 

Merge data

# create sequence of dates
date_sequence <- seq.Date(min(df_us_prev$date),
                          max(df_us_prev$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# join data frames 
df_us_prev <- df_us_prev %>%
  plyr::join(df_us_ctrl, by='county_fips') %>% 
  plyr::join(df_us_pers, by='county_fips') %>%
  merge(df_dates, by='date') %>% 
  arrange(county_fips, date)

df_us_prev

# create sequence of dates
date_sequence <- seq.Date(min(df_us_socdist$date),
                          as.Date('2020-04-28'), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# join data frames 
df_us_socdist <- df_us_socdist %>%
  plyr::join(df_us_ctrl, by='county_fips') %>% 
  plyr::join(df_us_pers, by='county_fips') %>%
  inner_join(df_dates, by='date') %>% 
  arrange(county_fips, date)

fips_complete <- df_us_socdist %>% 
  group_by(county_fips) %>% 
  summarize(n = n()) %>% 
  filter(! n<max(n)) %>% .$county_fips

df_us_socdist <- df_us_socdist %>%
  filter(county_fips %in% fips_complete)

df_us_socdist

Control for weekend effect

Plot prevalence over time

df_us_prev %>% sample_n(20000) %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_prev %>% 
  mutate(prev_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm=T), 
                                    quantile(.[[i]], 0.95, na.rm=T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") + 
  ggtitle(i)

print(gg)
}

Plot social distancing single tile visited

df_us_socdist %>% sample_n(10000) %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_socdist %>% 
  mutate(dist_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm = T), 
                                    quantile(.[[i]], 0.95, na.rm = T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(dist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~dist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_us_socdist %>% sample_n(10000) %>%
  ggplot(aes(x=time, y=socdist_tiles)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_socdist %>% 
  mutate(dist_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm = T), 
                                    quantile(.[[i]], 0.95, na.rm = T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(dist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_tiles)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~dist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_us_socdist <- df_us_socdist %>% 
  mutate(socdist_single_tile = socdist_single_tile_clean,
         socdist_tiles = socdist_tiles_clean) %>% 
  select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)

Correlations

df_us_prev %>% select(-time, -date) %>% 
  group_by(county_fips) %>%
  summarize_if(is.numeric, mean) %>% 
  select(-county_fips) %>%
  cor(use='pairwise.complete.obs') %>% 
  round(3) %>% as.data.frame()

df_us_socdist %>% select(-time, -date) %>% 
  group_by(county_fips) %>%
  summarize_if(is.numeric, mean) %>% 
  select(-county_fips) %>%
  cor(use='pairwise.complete.obs') %>% 
  round(3) %>% as.data.frame()

NA

Rescale Data

lvl2_scaled <- df_us_prev %>% 
  select(-time, -date, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_prev %>% select(county_fips, time, rate_day)

df_us_prev_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

lvl2_scaled <- df_us_socdist %>% 
  select(-time, -date, -socdist_tiles, -socdist_single_tile) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_socdist %>% 
  select(county_fips, time, socdist_single_tile, socdist_tiles)

df_us_socdist_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

Predict Prevalence

Extract first day of covid outbreak

# get onset day
df_us_onset_prev <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time)) %>%
  mutate(county_fips = as.character(county_fips))
  
# merge with county data
df_us_onset_prev <- df_us_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  mutate(county_fips = as.character(county_fips)) %>%
  left_join(df_us_onset_prev, by = 'county_fips')

# handle censored data
df_us_onset_prev <- df_us_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_us_prev$date)))+1))

Extract slopes

# cut time series before onset
df_us_prev_scaled <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_us_prev_scaled <- df_us_prev_scaled %>%
  group_by(county_fips) %>%
  filter(n() == 30) %>%
  ungroup()

# extract slope prevalence
df_us_slope_prev <- df_us_prev_scaled %>% split(.$county_fips) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('county_fips') %>% 
  rename(slope_prev = '.')
  
# merge with county data
df_us_slope_prev <- df_us_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_us_slope_prev, by = 'county_fips') %>%
  drop_na()

# handle censored data
df_us_slope_prev_cens <- df_us_slope_prev %>% 
  filter(slope_prev>=0) %>% 
  mutate(slope_prev = scale(slope_prev))

# standardize slopes
df_us_slope_prev <- df_us_slope_prev %>% 
  mutate(slope_prev = scale(slope_prev))

Explore distributions

df_us_onset_prev %>% ggplot(aes(onset_prev)) + geom_histogram(bins = 100)

df_us_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)

df_us_slope_prev_cens %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)

df_us_slope_prev %>% select(slope_prev) %>% summary()

    slope_prev.V1   
 Min.   :-0.437156  
 1st Qu.:-0.298897  
 Median :-0.216650  
 Mean   : 0.000000  
 3rd Qu.:-0.026712  
 Max.   :23.459218  

df_us_slope_prev_cens %>% select(slope_prev) %>% summary()

    slope_prev.V1   
 Min.   :-0.381258  
 1st Qu.:-0.300026  
 Median :-0.217861  
 Mean   : 0.000000  
 3rd Qu.:-0.029197  
 Max.   :23.334842  

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_us_onset_prev)
cox_onset_prev %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_us_onset_prev)

  n= 2469, number of events= 2429 
   (440 observations deleted due to missingness)

           coef exp(coef) se(coef)      z Pr(>|z|)    
pers_o  0.25801   1.29435  0.01953 13.214  < 2e-16 ***
pers_c  0.01242   1.01250  0.02440  0.509   0.6107    
pers_e  0.04044   1.04127  0.01961  2.062   0.0392 *  
pers_a  0.02804   1.02844  0.02404  1.167   0.2434    
pers_n -0.13832   0.87082  0.02214 -6.248 4.17e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.2944     0.7726    1.2458    1.3448
pers_c    1.0125     0.9877    0.9652    1.0621
pers_e    1.0413     0.9604    1.0020    1.0821
pers_a    1.0284     0.9723    0.9811    1.0781
pers_n    0.8708     1.1483    0.8338    0.9094

Concordance= 0.658  (se = 0.006 )
Likelihood ratio test= 280.5  on 5 df,   p=<2e-16
Wald test            = 282.1  on 5 df,   p=<2e-16
Score (logrank) test = 280.1  on 5 df,   p=<2e-16

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_distance + republican + medage + male + popdens +
                               manufact + tourism + academics + medinc + physician_pc,
                             data = df_us_onset_prev)
cox_onset_prev_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_distance + republican + medage + 
    male + popdens + manufact + tourism + academics + medinc + 
    physician_pc, data = df_us_onset_prev)

  n= 2464, number of events= 2424 
   (445 observations deleted due to missingness)

                      coef exp(coef)  se(coef)       z Pr(>|z|)    
pers_o            0.219404  1.245334  0.024013   9.137  < 2e-16 ***
pers_c            0.074751  1.077615  0.026097   2.864 0.004179 ** 
pers_e            0.016376  1.016510  0.020873   0.785 0.432729    
pers_a            0.042856  1.043787  0.026698   1.605 0.108444    
pers_n            0.001372  1.001373  0.025862   0.053 0.957701    
airport_distance -0.293113  0.745938  0.029875  -9.811  < 2e-16 ***
republican       -0.181505  0.834014  0.025588  -7.093 1.31e-12 ***
medage           -0.277857  0.757405  0.022256 -12.484  < 2e-16 ***
male             -0.145115  0.864923  0.028091  -5.166 2.39e-07 ***
popdens          -0.090498  0.913477  0.024988  -3.622 0.000293 ***
manufact          0.046112  1.047192  0.024014   1.920 0.054827 .  
tourism           0.068927  1.071358  0.021260   3.242 0.001186 ** 
academics         0.052582  1.053989  0.036863   1.426 0.153748    
medinc            0.308101  1.360838  0.032191   9.571  < 2e-16 ***
physician_pc     -0.107555  0.898027  0.022305  -4.822 1.42e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              1.2453     0.8030    1.1881    1.3053
pers_c              1.0776     0.9280    1.0239    1.1342
pers_e              1.0165     0.9838    0.9758    1.0590
pers_a              1.0438     0.9580    0.9906    1.0999
pers_n              1.0014     0.9986    0.9519    1.0534
airport_distance    0.7459     1.3406    0.7035    0.7909
republican          0.8340     1.1990    0.7932    0.8769
medage              0.7574     1.3203    0.7251    0.7912
male                0.8649     1.1562    0.8186    0.9139
popdens             0.9135     1.0947    0.8698    0.9593
manufact            1.0472     0.9549    0.9990    1.0977
tourism             1.0714     0.9334    1.0276    1.1169
academics           1.0540     0.9488    0.9805    1.1330
medinc              1.3608     0.7348    1.2776    1.4495
physician_pc        0.8980     1.1136    0.8596    0.9382

Concordance= 0.735  (se = 0.005 )
Likelihood ratio test= 998.1  on 15 df,   p=<2e-16
Wald test            = 1107  on 15 df,   p=<2e-16
Score (logrank) test = 1225  on 15 df,   p=<2e-16

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_us_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6944 -0.2924 -0.1891 -0.0187 23.3421 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.0005309  0.0209836  -0.025  0.97982   
pers_o      -0.0483148  0.0239424  -2.018  0.04371 * 
pers_c      -0.0198864  0.0290107  -0.685  0.49311   
pers_e      -0.0167807  0.0240723  -0.697  0.48581   
pers_a       0.0820008  0.0295572   2.774  0.00558 **
pers_n      -0.0468186  0.0271427  -1.725  0.08468 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9959 on 2254 degrees of freedom
Multiple R-squared:  0.01033,   Adjusted R-squared:  0.008138 
F-statistic: 4.707 on 5 and 2254 DF,  p-value: 0.0002796

lm_slope_prev %>% confint(level=0.9)

                    5 %         95 %
(Intercept) -0.03506004  0.033998144
pers_o      -0.08771282 -0.008916834
pers_c      -0.06762437  0.027851544
pers_e      -0.05639238  0.022830892
pers_a       0.03336354  0.130637980
pers_n      -0.09148282 -0.002154407

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc,
                         data = df_us_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc, 
    data = df_us_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7621 -0.3059 -0.1365  0.0537 22.0073 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -0.01144    0.02173  -0.526  0.59878    
pers_o           -0.01455    0.02781  -0.523  0.60097    
pers_c           -0.01846    0.02903  -0.636  0.52502    
pers_e           -0.01202    0.02408  -0.499  0.61775    
pers_a            0.08341    0.03049   2.735  0.00628 ** 
pers_n           -0.02391    0.02928  -0.817  0.41420    
airport_distance  0.04606    0.02877   1.601  0.10952    
republican       -0.10621    0.02598  -4.088 4.50e-05 ***
medage           -0.06991    0.02295  -3.047  0.00234 ** 
male              0.07055    0.02750   2.566  0.01036 *  
popdens           0.05905    0.02152   2.744  0.00611 ** 
manufact          0.08320    0.02525   3.295  0.00100 ** 
tourism           0.01052    0.02769   0.380  0.70401    
academics        -0.09219    0.03968  -2.323  0.02025 *  
medinc            0.09841    0.03352   2.936  0.00336 ** 
physician_pc      0.10935    0.02297   4.760 2.06e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9776 on 2244 degrees of freedom
Multiple R-squared:  0.05072,   Adjusted R-squared:  0.04437 
F-statistic: 7.993 on 15 and 2244 DF,  p-value: < 2.2e-16

lm_slope_prev_ctrl %>% confint(level=0.9)

                          5 %        95 %
(Intercept)      -0.047196737  0.02432453
pers_o           -0.060317938  0.03122009
pers_c           -0.066226814  0.02931540
pers_e           -0.051636468  0.02760241
pers_a            0.033234491  0.13358639
pers_n           -0.072089209  0.02426713
airport_distance -0.001281029  0.09340781
republican       -0.148955618 -0.06345988
medage           -0.107671822 -0.03215158
male              0.025301526  0.11579161
popdens           0.023640062  0.09445494
manufact          0.041644874  0.12474968
tourism          -0.035046903  0.05609068
academics        -0.157491478 -0.02689724
medinc            0.043247412  0.15356267
physician_pc      0.071549935  0.14715975

Predict prevalence slopes with censored regression

# predict slopes from personality
tobit_slope_prev <- censReg(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                            left = min(df_us_slope_prev_cens$slope_prev), right = Inf,
                         data = df_us_slope_prev_cens)
tobit_slope_prev %>% summary()

Call:
censReg(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, left = min(df_us_slope_prev_cens$slope_prev), right = Inf, 
    data = df_us_slope_prev_cens)

Observations:
         Total  Left-censored     Uncensored Right-censored 
          2233             14           2219              0 

Coefficients:
             Estimate Std. error t value Pr(> t)   
(Intercept) -0.005082   0.021189  -0.240 0.81046   
pers_o      -0.047752   0.024299  -1.965 0.04939 * 
pers_c      -0.017213   0.029440  -0.585 0.55875   
pers_e      -0.016517   0.024479  -0.675 0.49985   
pers_a       0.081401   0.030120   2.703 0.00688 **
pers_n      -0.047847   0.027592  -1.734 0.08290 . 
logSigma    -0.001513   0.015018  -0.101 0.91975   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Newton-Raphson maximisation, 2 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -3157.046 on 7 Df

tobit_slope_prev %>% confint(level=0.9)

                    5 %         95 %
(Intercept) -0.03993518  0.029771187
pers_o      -0.08771989 -0.007784003
pers_c      -0.06563724  0.031210354
pers_e      -0.05678204  0.023747982
pers_a       0.03185851  0.130943213
pers_n      -0.09323144 -0.002461866
logSigma    -0.02621535  0.023189436

# predict slopes from personality with controls
tobit_slope_prev_ctrl <- censReg(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc,
                           left = min(df_us_slope_prev_cens$slope_prev), right = Inf,
                         data = df_us_slope_prev_cens)
tobit_slope_prev_ctrl %>% summary()

Call:
censReg(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc, 
    left = min(df_us_slope_prev_cens$slope_prev), right = Inf, 
    data = df_us_slope_prev_cens)

Observations:
         Total  Left-censored     Uncensored Right-censored 
          2233             14           2219              0 

Coefficients:
                 Estimate Std. error t value  Pr(> t)    
(Intercept)      -0.01735    0.02201  -0.788 0.430514    
pers_o           -0.01439    0.02824  -0.509 0.610412    
pers_c           -0.01520    0.02940  -0.517 0.605235    
pers_e           -0.01174    0.02446  -0.480 0.631236    
pers_a            0.08337    0.03098   2.691 0.007128 ** 
pers_n           -0.02222    0.02974  -0.747 0.454893    
airport_distance  0.04792    0.02913   1.645 0.099979 .  
republican       -0.10787    0.02622  -4.114 3.88e-05 ***
medage           -0.07138    0.02327  -3.068 0.002158 ** 
male              0.06885    0.02761   2.494 0.012643 *  
popdens           0.06047    0.02166   2.792 0.005243 ** 
manufact          0.08489    0.02540   3.342 0.000833 ***
tourism           0.01388    0.02807   0.494 0.620985    
academics        -0.09198    0.03987  -2.307 0.021045 *  
medinc            0.09873    0.03369   2.931 0.003380 ** 
physician_pc      0.11114    0.02330   4.769 1.85e-06 ***
logSigma         -0.02252    0.01502  -1.499 0.133768    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Newton-Raphson maximisation, 3 iterations
Return code 1: gradient close to zero
Log-likelihood: -3110.339 on 17 Df

tobit_slope_prev_ctrl %>% confint(level=0.9)

                           5 %         95 %
(Intercept)      -5.354314e-02  0.018848670
pers_o           -6.084617e-02  0.032066833
pers_c           -6.355657e-02  0.033162950
pers_e           -5.196628e-02  0.028488571
pers_a            3.240711e-02  0.134335675
pers_n           -7.113228e-02  0.026690043
airport_distance  2.991424e-06  0.095841919
republican       -1.509982e-01 -0.064747430
medage           -1.096460e-01 -0.033104079
male              2.343624e-02  0.114267420
popdens           2.484103e-02  0.096096473
manufact          4.310256e-02  0.126669193
tourism          -3.228815e-02  0.060043758
academics        -1.575597e-01 -0.026405710
medinc            4.331968e-02  0.154134078
physician_pc      7.280430e-02  0.149466568
logSigma         -4.721970e-02  0.002184288

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp

          pers_o           pers_c           pers_e           pers_a 
    0.0069588100     0.0053233415     0.0009241402     0.0336522004 
          pers_n airport_distance       republican           medage 
    0.0101345184     0.0088034043     0.0113297404     0.0409948903 
            male          popdens         manufact          tourism 
    0.0159517791     0.0237333099     0.0559568844     0.0099654006 
       academics           medinc     physician_pc 
    0.0160422150     0.0067803089    -0.0006294380 

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond

          pers_o           pers_c           pers_e           pers_a 
    0.0030020160     0.0186830072     0.0003241834     0.0312488347 
          pers_n airport_distance       republican           medage 
    0.0131618928     0.0212862058     0.0088770492     0.0536994750 
            male          popdens         manufact          tourism 
    0.0148531988     0.0239457182     0.0609858257     0.0094635969 
       academics           medinc     physician_pc 
    0.0219869043     0.0204296385     0.0013880483 

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
fips_complete <- df_us_socdist_scaled %>% 
  group_by(county_fips) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$county_fips

# run changepoint analysis
df_us_socdist_cpt_results <- df_us_socdist_scaled %>% 
  select(county_fips, socdist_single_tile) %>%
  filter(county_fips %in% fips_complete) %>% 
  split(.$county_fips) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

df_us_socdist_cpt_results_2 <- df_us_socdist_scaled %>% 
  select(county_fips, socdist_tiles) %>%
  filter(county_fips %in% fips_complete) %>% 
  split(.$county_fips) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_tiles),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_us_socdist_cpt_day <- df_us_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('county_fips')

df_us_socdist_cpt_day_2 <- df_us_socdist_cpt_results_2 %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist_2 = '.') %>%
  rownames_to_column('county_fips')

# calculate mean differences
df_us_socdist_cpt_mean_diff <- df_us_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('county_fips')

df_us_socdist_cpt_mean_diff_2 <- df_us_socdist_cpt_results_2 %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[1]-.[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist_2 = '.') %>%
  rownames_to_column('county_fips')

# merge with county data
df_us_cpt_socdist <- df_us_socdist_scaled %>% 
  select(-time, -socdist_single_tile, -socdist_tiles) %>%
  distinct() %>% 
  left_join(df_us_socdist_cpt_day, by='county_fips') %>%
  left_join(df_us_socdist_cpt_day_2, by='county_fips') %>%
  left_join(df_us_socdist_cpt_mean_diff, by='county_fips') %>%
  left_join(df_us_socdist_cpt_mean_diff_2, by='county_fips') %>%
  left_join(select(df_us_onset_prev, county_fips, onset_prev), by='county_fips') %>%
  left_join(select(df_us_slope_prev, county_fips, slope_prev), by='county_fips') 

# standardize mean/var differences
df_us_cpt_socdist <- df_us_cpt_socdist %>% 
  mutate(mean_diff_socdist = scale(mean_diff_socdist),
         mean_diff_socdist_2 = scale(mean_diff_socdist_2))

# handle censored data
df_us_cpt_socdist <- df_us_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), as.numeric(diff(range(df_us_socdist$date))), cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= as.numeric(diff(range(df_us_socdist$date))), 0, 1))

df_us_cpt_socdist$cpt_day_socdist %>% hist()

df_us_cpt_socdist$mean_diff_socdist %>% hist()

df_us_cpt_socdist$cpt_day_socdist_2 %>% hist()

df_us_cpt_socdist$mean_diff_socdist_2 %>% hist()

for(i in head(df_us_socdist_cpt_results, 5)){
  plot(i)
}

NA

for(i in head(df_us_socdist_cpt_results_2, 5)){
  plot(i)
}

NA

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_us_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_us_cpt_socdist)

  n= 2380, number of events= 2379 

            coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o -0.189086  0.827715  0.023503 -8.045 8.61e-16 ***
pers_c -0.033629  0.966931  0.027365 -1.229  0.21912    
pers_e  0.004178  1.004187  0.023905  0.175  0.86126    
pers_a -0.084190  0.919256  0.028674 -2.936  0.00332 ** 
pers_n  0.060781  1.062667  0.027280  2.228  0.02588 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    0.8277     1.2081    0.7905    0.8667
pers_c    0.9669     1.0342    0.9164    1.0202
pers_e    1.0042     0.9958    0.9582    1.0524
pers_a    0.9193     1.0878    0.8690    0.9724
pers_n    1.0627     0.9410    1.0073    1.1210

Concordance= 0.593  (se = 0.007 )
Likelihood ratio test= 121.9  on 5 df,   p=<2e-16
Wald test            = 124.3  on 5 df,   p=<2e-16
Score (logrank) test = 123.9  on 5 df,   p=<2e-16

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                           onset_prev + slope_prev ,
                  data = df_us_cpt_socdist)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_distance + republican + 
    medage + male + popdens + manufact + tourism + academics + 
    medinc + physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

  n= 2196, number of events= 2195 
   (184 observations deleted due to missingness)

                       coef  exp(coef)   se(coef)      z Pr(>|z|)    
pers_o           -0.0923731  0.9117649  0.0296366 -3.117 0.001828 ** 
pers_c           -0.0557289  0.9457955  0.0290405 -1.919 0.054983 .  
pers_e            0.0024032  1.0024061  0.0260018  0.092 0.926360    
pers_a           -0.0834917  0.9198987  0.0321849 -2.594 0.009483 ** 
pers_n            0.0218463  1.0220867  0.0318519  0.686 0.492794    
airport_distance  0.0323170  1.0328449  0.0270192  1.196 0.231668    
republican        0.0274670  1.0278477  0.0270193  1.017 0.309358    
medage            0.0168080  1.0169500  0.0221423  0.759 0.447799    
male              0.0123963  1.0124734  0.0233575  0.531 0.595614    
popdens          -0.2671107  0.7655883  0.0721589 -3.702 0.000214 ***
manufact          0.0983611  1.1033612  0.0242468  4.057 4.98e-05 ***
tourism          -0.0858071  0.9177712  0.0265987 -3.226 0.001255 ** 
academics         0.1451860  1.1562546  0.0431797  3.362 0.000773 ***
medinc           -0.0752502  0.9275114  0.0360780 -2.086 0.037000 *  
physician_pc      0.0001824  1.0001824  0.0253754  0.007 0.994264    
onset_prev        0.0077552  1.0077853  0.0032047  2.420 0.015522 *  
slope_prev       -0.1210372  0.8860010  0.0317937 -3.807 0.000141 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              0.9118     1.0968    0.8603    0.9663
pers_c              0.9458     1.0573    0.8935    1.0012
pers_e              1.0024     0.9976    0.9526    1.0548
pers_a              0.9199     1.0871    0.8637    0.9798
pers_n              1.0221     0.9784    0.9602    1.0879
airport_distance    1.0328     0.9682    0.9796    1.0890
republican          1.0278     0.9729    0.9748    1.0837
medage              1.0170     0.9833    0.9738    1.0621
male                1.0125     0.9877    0.9672    1.0599
popdens             0.7656     1.3062    0.6646    0.8819
manufact            1.1034     0.9063    1.0522    1.1571
tourism             0.9178     1.0896    0.8712    0.9669
academics           1.1563     0.8649    1.0624    1.2584
medinc              0.9275     1.0782    0.8642    0.9955
physician_pc        1.0002     0.9998    0.9517    1.0512
onset_prev          1.0078     0.9923    1.0015    1.0141
slope_prev          0.8860     1.1287    0.8325    0.9430

Concordance= 0.615  (se = 0.008 )
Likelihood ratio test= 203.7  on 17 df,   p=<2e-16
Wald test            = 184.4  on 17 df,   p=<2e-16
Score (logrank) test = 188.3  on 17 df,   p=<2e-16

Linear models predicting mean differences

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_us_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7658 -0.5742 -0.0819  0.4870  6.0744 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.168e-16  1.808e-02   0.000        1    
pers_o       3.448e-01  1.961e-02  17.583  < 2e-16 ***
pers_c      -1.622e-01  2.453e-02  -6.610 4.72e-11 ***
pers_e       1.321e-01  2.017e-02   6.549 7.08e-11 ***
pers_a      -1.376e-01  2.505e-02  -5.492 4.39e-08 ***
pers_n      -1.104e-01  2.288e-02  -4.825 1.49e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8822 on 2374 degrees of freedom
Multiple R-squared:  0.2233,    Adjusted R-squared:  0.2217 
F-statistic: 136.5 on 5 and 2374 DF,  p-value: < 2.2e-16

lm_meandiff_socdist %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.02975669  0.02975669
pers_o       0.31251542  0.37704681
pers_c      -0.20255000 -0.12180865
pers_e       0.09890660  0.16529010
pers_a      -0.17877385 -0.09634921
pers_n      -0.14807707 -0.07276741

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n + 
                              airport_distance + republican + medage + male + popdens + 
                              manufact + tourism + academics + medinc + physician_pc +
                              onset_prev + slope_prev ,
                            data = df_us_cpt_socdist)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc + 
    onset_prev + slope_prev, data = df_us_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.18454 -0.37354 -0.04018  0.32901  2.71979 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.933063   0.119257   7.824 7.92e-15 ***
pers_o            0.050890   0.016370   3.109 0.001903 ** 
pers_c           -0.055675   0.017167  -3.243 0.001200 ** 
pers_e            0.059599   0.014328   4.160 3.31e-05 ***
pers_a           -0.042432   0.018380  -2.309 0.021061 *  
pers_n            0.131964   0.017479   7.550 6.36e-14 ***
airport_distance -0.039999   0.013904  -2.877 0.004057 ** 
republican       -0.260446   0.015120 -17.225  < 2e-16 ***
medage            0.099828   0.013375   7.464 1.21e-13 ***
male              0.043997   0.013892   3.167 0.001561 ** 
popdens           0.234350   0.031231   7.504 8.99e-14 ***
manufact          0.009840   0.014850   0.663 0.507633    
tourism           0.090837   0.014939   6.081 1.41e-09 ***
academics         0.161546   0.024473   6.601 5.12e-11 ***
medinc            0.378549   0.020675  18.309  < 2e-16 ***
physician_pc     -0.004576   0.013600  -0.336 0.736576    
onset_prev       -0.013586   0.001764  -7.703 2.00e-14 ***
slope_prev        0.042580   0.012512   3.403 0.000678 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5709 on 2178 degrees of freedom
  (184 observations deleted due to missingness)
Multiple R-squared:  0.6642,    Adjusted R-squared:  0.6616 
F-statistic: 253.5 on 17 and 2178 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl %>% confint(level=0.9)

                         5 %        95 %
(Intercept)       0.73681930  1.12930762
pers_o            0.02395230  0.07782772
pers_c           -0.08392327 -0.02742598
pers_e            0.03602151  0.08317608
pers_a           -0.07267787 -0.01218633
pers_n            0.10320155  0.16072607
airport_distance -0.06287822 -0.01711884
republican       -0.28532668 -0.23556552
medage            0.07781790  0.12183798
male              0.02113686  0.06685733
popdens           0.18295701  0.28574228
manufact         -0.01459616  0.03427626
tourism           0.06625519  0.11541954
academics         0.12127362  0.20181810
medinc            0.34452716  0.41257115
physician_pc     -0.02695577  0.01780451
onset_prev       -0.01648877 -0.01068414
slope_prev        0.02199129  0.06316910

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                                  airport_distance + republican + medage + male + popdens + 
                                  manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_cpt_socdist,
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

          pers_o           pers_c           pers_e           pers_a 
    1.614037e-02     2.423785e-02     2.732908e-05     3.101524e-02 
          pers_n airport_distance       republican           medage 
    9.194875e-03     6.062698e-03     1.533444e-01     1.409992e-02 
            male          popdens         manufact          tourism 
    6.976691e-03     1.440436e-01     7.407710e-03     1.852098e-02 
       academics           medinc     physician_pc 
    1.529007e-01     2.018907e-01     2.242878e-03 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

          pers_o           pers_c           pers_e           pers_a 
    0.0173120419     0.0268238558     0.0007997082     0.0325974638 
          pers_n airport_distance       republican           medage 
    0.0127582857     0.0065079941     0.1673575442     0.0147002262 
            male          popdens         manufact          tourism 
    0.0066593632     0.1439615291     0.0075370064     0.0190248211 
       academics           medinc     physician_pc 
    0.1686254063     0.2106451498     0.0026035500 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist_2 <- coxph(Surv(cpt_day_socdist_2, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_us_cpt_socdist)
cox_cpt_socdist_2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist_2, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_us_cpt_socdist)

  n= 2380, number of events= 2379 

           coef exp(coef) se(coef)      z Pr(>|z|)    
pers_o -0.14730   0.86303  0.02283 -6.453 1.10e-10 ***
pers_c -0.11744   0.88919  0.02688 -4.369 1.25e-05 ***
pers_e  0.03018   1.03064  0.02355  1.282  0.19998    
pers_a  0.03223   1.03276  0.02720  1.185  0.23607    
pers_n  0.07865   1.08183  0.02638  2.982  0.00287 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    0.8630     1.1587    0.8253    0.9025
pers_c    0.8892     1.1246    0.8436    0.9373
pers_e    1.0306     0.9703    0.9841    1.0793
pers_a    1.0328     0.9683    0.9791    1.0893
pers_n    1.0818     0.9244    1.0273    1.1392

Concordance= 0.598  (se = 0.008 )
Likelihood ratio test= 109.1  on 5 df,   p=<2e-16
Wald test            = 112.2  on 5 df,   p=<2e-16
Score (logrank) test = 111.5  on 5 df,   p=<2e-16

# predict hazard from personality with controls
cox_cpt_socdist_ctrl_2 <- coxph(Surv(cpt_day_socdist_2, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                           onset_prev + slope_prev ,
                  data = df_us_cpt_socdist)
cox_cpt_socdist_ctrl_2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist_2, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_distance + republican + 
    medage + male + popdens + manufact + tourism + academics + 
    medinc + physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

  n= 2196, number of events= 2195 
   (184 observations deleted due to missingness)

                      coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o           -0.103657  0.901534  0.029764 -3.483 0.000497 ***
pers_c           -0.120681  0.886317  0.029377 -4.108 3.99e-05 ***
pers_e            0.035021  1.035642  0.026256  1.334 0.182259    
pers_a            0.014626  1.014733  0.031050  0.471 0.637614    
pers_n            0.046338  1.047428  0.031442  1.474 0.140544    
airport_distance  0.005353  1.005367  0.025734  0.208 0.835217    
republican       -0.006681  0.993342  0.026860 -0.249 0.803578    
medage            0.010936  1.010996  0.022736  0.481 0.630517    
male             -0.010455  0.989599  0.023257 -0.450 0.653034    
popdens          -0.233974  0.791383  0.066998 -3.492 0.000479 ***
manufact          0.154049  1.166548  0.025299  6.089 1.14e-09 ***
tourism          -0.068794  0.933519  0.026663 -2.580 0.009876 ** 
academics         0.128110  1.136678  0.041308  3.101 0.001926 ** 
medinc           -0.011745  0.988323  0.034107 -0.344 0.730568    
physician_pc     -0.030768  0.969701  0.024904 -1.235 0.216659    
onset_prev        0.004743  1.004754  0.003206  1.479 0.139121    
slope_prev       -0.102630  0.902461  0.030615 -3.352 0.000801 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              0.9015     1.1092    0.8504    0.9557
pers_c              0.8863     1.1283    0.8367    0.9388
pers_e              1.0356     0.9656    0.9837    1.0903
pers_a              1.0147     0.9855    0.9548    1.0784
pers_n              1.0474     0.9547    0.9848    1.1140
airport_distance    1.0054     0.9947    0.9559    1.0574
republican          0.9933     1.0067    0.9424    1.0470
medage              1.0110     0.9891    0.9669    1.0571
male                0.9896     1.0105    0.9455    1.0358
popdens             0.7914     1.2636    0.6940    0.9024
manufact            1.1665     0.8572    1.1101    1.2259
tourism             0.9335     1.0712    0.8860    0.9836
academics           1.1367     0.8798    1.0483    1.2325
medinc              0.9883     1.0118    0.9244    1.0566
physician_pc        0.9697     1.0312    0.9235    1.0182
onset_prev          1.0048     0.9953    0.9985    1.0111
slope_prev          0.9025     1.1081    0.8499    0.9583

Concordance= 0.63  (se = 0.008 )
Likelihood ratio test= 197.1  on 17 df,   p=<2e-16
Wald test            = 188.2  on 17 df,   p=<2e-16
Score (logrank) test = 190.4  on 17 df,   p=<2e-16

Linear models predicting mean differences

lm_meandiff_socdist_2 <- lm(mean_diff_socdist_2 ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_cpt_socdist)
lm_meandiff_socdist_2 %>% summary()

Call:
lm(formula = mean_diff_socdist_2 ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_us_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9433  -0.3894   0.0131   0.4666   3.5265 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.037e-16  1.935e-02   0.000 1.000000    
pers_o       3.012e-01  2.098e-02  14.355  < 2e-16 ***
pers_c      -9.312e-02  2.625e-02  -3.547 0.000397 ***
pers_e       6.867e-02  2.158e-02   3.181 0.001484 ** 
pers_a      -2.872e-02  2.680e-02  -1.072 0.284018    
pers_n      -2.435e-02  2.449e-02  -0.994 0.320129    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.944 on 2374 degrees of freedom
Multiple R-squared:  0.1107,    Adjusted R-squared:  0.1089 
F-statistic: 59.13 on 5 and 2374 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_2 %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.03184033  0.03184033
pers_o       0.26667433  0.33572435
pers_c      -0.13632049 -0.04992544
pers_e       0.03315103  0.10418286
pers_a      -0.07281565  0.01538056
pers_n      -0.06464033  0.01594268

lm_meandiff_socdist_ctrl_2 <- lm(mean_diff_socdist_2 ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n + 
                              airport_distance + republican + medage + male + popdens + 
                              manufact + tourism + academics + medinc + physician_pc +
                              onset_prev + slope_prev ,
                            data = df_us_cpt_socdist)
lm_meandiff_socdist_ctrl_2 %>% summary()

Call:
lm(formula = mean_diff_socdist_2 ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_distance + republican + medage + 
    male + popdens + manufact + tourism + academics + medinc + 
    physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.7181  -0.2980   0.0295   0.3613   3.5496 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.69264    0.17509   3.956 7.87e-05 ***
pers_o            0.10029    0.02403   4.173 3.13e-05 ***
pers_c           -0.03809    0.02520  -1.511  0.13090    
pers_e            0.02990    0.02104   1.421  0.15543    
pers_a            0.03920    0.02699   1.453  0.14647    
pers_n            0.13351    0.02566   5.203 2.15e-07 ***
airport_distance -0.11922    0.02041  -5.840 6.00e-09 ***
republican       -0.16546    0.02220  -7.453 1.31e-13 ***
medage            0.05129    0.01964   2.612  0.00906 ** 
male              0.01352    0.02040   0.663  0.50753    
popdens           0.05980    0.04585   1.304  0.19229    
manufact         -0.10282    0.02180  -4.716 2.56e-06 ***
tourism           0.09773    0.02193   4.456 8.77e-06 ***
academics        -0.03454    0.03593  -0.961  0.33649    
medinc            0.30234    0.03036   9.960  < 2e-16 ***
physician_pc     -0.03573    0.01997  -1.789  0.07371 .  
onset_prev       -0.01011    0.00259  -3.905 9.70e-05 ***
slope_prev        0.03842    0.01837   2.091  0.03662 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8382 on 2178 degrees of freedom
  (184 observations deleted due to missingness)
Multiple R-squared:  0.2922,    Adjusted R-squared:  0.2867 
F-statistic:  52.9 on 17 and 2178 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl_2 %>% confint(level=0.9)

                          5 %         95 %
(Intercept)       0.404509277  0.980759910
pers_o            0.060739497  0.139839292
pers_c           -0.079561227  0.003387984
pers_e           -0.004721244  0.064511015
pers_a           -0.005205988  0.083607586
pers_n            0.091283064  0.175740460
airport_distance -0.152808310 -0.085624474
republican       -0.201985450 -0.128926199
medage            0.018978426  0.083608629
male             -0.020044734  0.047081972
popdens          -0.015651171  0.135257974
manufact         -0.138694111 -0.066939714
tourism           0.061641628  0.133824640
academics        -0.093670255  0.024584991
medinc            0.252393233  0.352295292
physician_pc     -0.068586967 -0.002870009
onset_prev       -0.014374120 -0.005851771
slope_prev        0.008188419  0.068645614

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist_2 <- cforest(mean_diff_socdist_2 ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                                  airport_distance + republican + medage + male + popdens + 
                                  manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_cpt_socdist,
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist_2, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist_2, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

          pers_o           pers_c           pers_e           pers_a 
     0.039260552      0.012302722      0.002330488      0.011952030 
          pers_n airport_distance       republican           medage 
     0.006029187      0.020302199      0.080534580      0.009563549 
            male          popdens         manufact          tourism 
     0.005388067      0.108579172      0.036924193      0.016906413 
       academics           medinc     physician_pc 
     0.055033273      0.090556592      0.002284354 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

          pers_o           pers_c           pers_e           pers_a 
     0.043664416      0.011665283      0.001743354      0.014833439 
          pers_n airport_distance       republican           medage 
     0.008701287      0.020987417      0.081218508      0.012618969 
            male          popdens         manufact          tourism 
     0.006820099      0.110563441      0.036236182      0.016517248 
       academics           medinc     physician_pc 
     0.072596700      0.095937274      0.005176225 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Export data

us_list_results <- list(cox_onset_prev, cox_onset_prev_ctrl, 
     lm_slope_prev, lm_slope_prev_ctrl, 
     cox_cpt_socdist, cox_cpt_socdist_ctrl,
     lm_meandiff_socdist, lm_meandiff_socdist_ctrl,
     cox_cpt_socdist_2, cox_cpt_socdist_ctrl_2,
     lm_meandiff_socdist_2, lm_meandiff_socdist_ctrl_2)

results_names <- list('cox_onset_prev', 'cox_onset_prev_ctrl', 
     'lm_slope_prev', 'lm_slope_prev_ctrl', 
     'cox_cpt_socdist', 'cox_cpt_socdist_ctrl', 
     'lm_meandiff_socdist', 'lm_meandiff_socdist_ctrl',
     'cox_cpt_socdist_2', 'cox_cpt_socdist_ctrl_2',
     'lm_meandiff_socdist_2', 'lm_meandiff_socdist_ctrl_2')

names(us_list_results) <- results_names

save(us_list_results, file="us_list_results.RData")

write_csv(df_us_slope_prev, 'df_us_slope_prev.csv')
write_csv(df_us_cpt_socdist, 'df_us_cpt_socdist.csv')